wikiversitywikiaorg_ru-20200214-history
Программирование и научные вычисления на языке Python/§19
Хотя численные методы, некоторые из которых описаны в предыдущих уроках, и находят широкое применение, большей точностью и наглядностью обладают символьные вычисления, которые работают с математическими выражениями, собственно, как и предполагает математика, как с последовательностями символов. Системы, занимающиеся символьными вычислениями, называют также системами компьютерной алгебры. Примерами таких систем служат известные математические среды Maple, Mathcad, Mathematica, Maxima и т. д. В качестве инструмента символьных вычислений мы рассмотрим библиотеку SymPy. SymPy представляет собой библиотеку символьных вычислений, которая в конечной цели стремится стать полноценной системой компьютерной алгебры, сохраняя при этом как можно более простой код, ясный для понимания и дальнейших изменений и дополнений. SymPy написан исключительно на Python и не требует никаких других библиотек. Установка SymPy происходит аналогично другим продуктам, рассмотренным ранее. isympy Если вы уже используете для своих научных экспериментов IPython, то для вас окажется приятным тот факт, что в SymPy для него имеется обертка — isympy, запуск которой в консоли позволяет упростить некоторые стандартные действия относительно переменных и представлять данные на выходе в удобном для понимания формате записи: leo@leo:~$ isympy Python 2.6.5 console for SymPy 0.6.6 These commands were executed: >>> from __future__ import division >>> from sympy import * >>> x, y, z = symbols('xyz') >>> k, m, n = symbols('kmn', integer=True) >>> f, g, h = map(Function, 'fgh') Documentation can be found at http://sympy.org/ In 1: (1/cos(x)).series(x, 0, 10) Out1: 2 4 6 8 x 5⋅x 61⋅x 277⋅x 1 + ── + ──── + ───── + ────── + O(x**10) 2 24 720 8064 Далее в уроке, в зависимости от удобства представления, одни данные будут показаны как при работе в стандартном IDLE, другие — как при работе в isympy. Калькулятор SymPy содержит три встроенных числовых типа: Real, Rational и Integer. Класс Rational отвечает за дроби, любую пару целых чисел он определяет как числитель и знаменатель, далее над дробями можно производить типичные действия: умножать, делить, складывать и так далее, при этом например: >>> from sympy import * >>> a = Rational(1,2) >>> a 1/2 >>> a*2 1 >>> Rational(2)**50/Rational(10)**50 1/88817841970012523233890533447265625 # при этом: >>> 1/2 0 >>> 1.0/2 0.5 # можно использовать математические константы >>> pi**2 pi**2 >>> pi.evalf() 3.14159265358979 >>> (pi+exp(1)).evalf() 5.85987448204884 # а также и символ для бесконечности >>> oo > 99999 True >>> oo + 1 oo Символы В отличие от других систем компьютерной алгебры, в SymPy необходимо предварительно обозначить, какие переменные являются символьными: >>> from sympy import * >>> x = Symbol('x') >>> y = Symbol('y') Теперь вы можете играть с ними: >>> x+y+x-y 2*x >>> (x+y)**2 (x + y)**2 >>> ((x+y)**2).expand() 2*x*y + x**2 + y**2 И заменять их другими символами или числами, используя конструкцию subs(old, new): >>> ((x+y)**2).subs(x, 1) (1 + y)**2 >>> ((x+y)**2).subs(x, y) 4*y**2 Разложение на элементарные Для разложения дроби на элементарные используется apart(expr, x): In 1: 1/( (x+2)*(x+1) ) Out1: 1 ─────────────── (2 + x)*(1 + x) In 2: apart(1/( (x+2)*(x+1) ), x) Out2: 1 1 ───── - ───── 1 + x 2 + x In 3: (x+1)/(x-1) Out3: -(1 + x) ──────── 1 - x In 4: apart((x+1)/(x-1), x) Out4: 2 1 - ───── 1 - x Чтобы снова представить их в виде одной дроби together(expr, x): In 7: together(1/x + 1/y + 1/z) Out7: x*y + x*z + y*z ─────────────── x*y*z In 8: together(apart((x+1)/(x-1), x), x) Out8: -1 - x ────── 1 - x In 9: together(apart(1/( (x+2)*(x+1) ), x), x) Out9: 1 ─────────────── (2 + x)*(1 + x) Пределы Пределы в SymPy определяются очень просто. Чтобы найти предел функции f(x) при x -> 0, необходимо так и записать limit(f(x), x, 0): >>> from sympy import * >>> x=Symbol("x") >>> limit(sin(x)/x, x, 0) 1 Естественно, можно использовать и предел в бесконечности: >>> limit(x, x, oo) oo >>> limit(1/x, x, oo) 0 >>> limit(x**x, x, 0) 1 Дифференцирование Любое SymPy выражение можно продифференцировать, используя diff(func, var): >>> from sympy import * >>> x = Symbol('x') >>> diff(sin(x), x) cos(x) >>> diff(sin(2*x), x) 2*cos(2*x) >>> diff(tan(x), x) 1 + tan(x)**2 По определению производной вы можете проверить правильность решения: >>> limit((tan(x+y)-tan(x))/y, y, 0) 1 + tan(x)**2 Вторые, третьи и так далее производные можно найти с помощью diff(func, var, n): >>> diff(sin(2*x), x, 1) 2*cos(2*x) >>> diff(sin(2*x), x, 2) -4*sin(2*x) >>> diff(sin(2*x), x, 3) -8*cos(2*x) Разложение в ряд Используем метод .series(var, point, order): >>> from sympy import * >>> x = Symbol('x') >>> cos(x).series(x, 0, 10) 1 - x**2/2 + x**4/24 - x**6/720 + x**8/40320 + O(x**10) >>> (1/cos(x)).series(x, 0, 10) 1 + x**2/2 + 5*x**4/24 + 61*x**6/720 + 277*x**8/8064 + O(x**10) Другой простой пример: from sympy import Symbol, pprint x = Symbol("x") y = Symbol("y") e = 1/(x + y) s = e.series(x, 0, 5) print(s) pprint(s) дает нам после запуска: 1/y + x**2*y**(-3) + x**4*y**(-5) - x*y**(-2) - x**3*y**(-4) + O(x**5) 2 4 3 1 x x x x ─ + ── + ── - ── - ── + O(x**5) y 3 5 2 4 y y y y Интегрирование SymPy также позволяет брать неопределенные и определенные интегралы как от трансцендентных элементарных, так и от специальных функций: >>> from sympy import * >>> x, y = symbols('xy') Интегралы от элементарных функций: >>> integrate(6*x**5, x) x**6 >>> integrate(sin(x), x) -cos(x) >>> integrate(log(x), x) -x + x*log(x) >>> integrate(2*x + sinh(x), x) cosh(x) + x**2 И от специальных: >>> integrate(exp(-x**2)*erf(x), x) pi**(1/2)*erf(x)**2/4 Можно найти и определенный интеграл: >>> integrate(x**3, (x, -1, 1)) 0 >>> integrate(sin(x), (x, 0, pi/2)) 1 >>> integrate(cos(x), (x, -pi/2, pi/2)) 2 Или несобственный интеграл: >>> integrate(exp(-x), (x, 0, oo)) 1 >>> integrate(log(x), (x, 0, 1)) -1 Комплексные числа >>> from sympy import Symbol, exp, I >>> x = Symbol("x") >>> exp(I*x).expand() exp(I*x) >>> exp(I*x).expand(complex=True) I*exp(-im(x))*sin(re(x)) + cos(re(x))*exp(-im(x)) >>> x = Symbol("x", real=True) >>> exp(I*x).expand(complex=True) I*sin(x) + cos(x) Функции Тригонометрические In 1: sin(x+y).expand(trig=True) Out1: cos(x)*sin(y) + cos(y)*sin(x) In 2: cos(x+y).expand(trig=True) Out2: cos(x)*cos(y) - sin(x)*sin(y) In 3: sin(I*x) Out3: ⅈ*sinh(x) In 4: sinh(I*x) Out4: ⅈ*sin(x) In 5: asinh(I) Out5: π*ⅈ ─── 2 In 6: asinh(I*x) Out6: ⅈ*asin(x) In 15: sin(x).series(x, 0, 10) Out15: 3 5 7 9 x x x x x - ── + ─── - ──── + ────── + O(x**10) 6 120 5040 362880 In 16: sinh(x).series(x, 0, 10) Out16: 3 5 7 9 x x x x x + ── + ─── + ──── + ────── + O(x**10) 6 120 5040 362880 In 17: asin(x).series(x, 0, 10) Out17: 3 5 7 9 x 3*x 5*x 35*x x + ── + ──── + ──── + ───── + O(x**10) 6 40 112 1152 In 18: asinh(x).series(x, 0, 10) Out18: 3 5 7 9 x 3*x 5*x 35*x x - ── + ──── - ──── + ───── + O(x**10) 6 40 112 1152 Сферические гармоники In 1: from sympy.abc import theta, phi In 2: Ylm(1, 0, theta, phi) Out2: ⎽⎽⎽ ╲╱ 3 *cos(θ) ──────────── ⎽⎽⎽ 2*╲╱ π In 3: Ylm(1, 1, theta, phi) Out3: ⎽⎽⎽ ⅈ*φ -╲╱ 6 *│sin(θ)│*ℯ ──────────────────── ⎽⎽⎽ 4*╲╱ π In 4: Ylm(2, 1, theta, phi) Out4: ⎽⎽⎽⎽ ⅈ*φ -╲╱ 30 *│sin(θ)│*cos(θ)*ℯ ──────────────────────────── ⎽⎽⎽ 4*╲╱ π Факториалы и гамма-функции In 1: x = Symbol("x") In 2: y = Symbol("y", integer=True) In 3: factorial(x) Out3: Γ(1 + x) In 4: factorial(y) Out4: y! In 5: factorial(x).series(x, 0, 3) Out5: 2 2 2 2 x *EulerGamma π *x 1 - x*EulerGamma + ────────────── + ───── + O(x**3) 2 12 Дзета-функция Римана In 18: zeta(4, x) Out18: ζ(4, x) In 19: zeta(4, 1) Out19: 4 π ── 90 In 20: zeta(4, 2) Out20: 4 π -1 + ── 90 In 21: zeta(4, 3) Out21: 4 17 π - ── + ── 16 90 Многочлены Полиномы Чебышева, Лежандра, Эрмита: In 1: chebyshevt(2, x) Out1: 2 -1 + 2*x In 2: chebyshevt(4, x) Out2: 2 4 1 - 8*x + 8*x In 3: legendre(2, x) Out3: 2 3*x -1/2 + ──── 2 In 4: legendre(8, x) Out4: 2 4 6 8 35 315*x 3465*x 3003*x 6435*x ─── - ────── + ─────── - ─────── + ─────── 128 32 64 32 128 In 5: assoc_legendre(2, 1, x) Out5: ⎽⎽⎽⎽⎽⎽⎽⎽ ╱ 2 -3*x*╲╱ 1 - x In 6: assoc_legendre(2, 2, x) Out6: 2 3 - 3*x In 7: hermite(3, x) Out7: 3 -12*x + 8*x Дифференциальные уравнения In 1: f(x).diff(x, x) + f(x) Out1: 2 d ─────(f(x)) + f(x) dx dx In 2: dsolve(f(x).diff(x, x) + f(x), f(x)) Out2: C₁*sin(x) + C₂*cos(x) Алгебраические уравнения In 3: solve(x**4 - 1, x) Out3: 1, -1, -ⅈ In 4: solve(+ 5*y - 2, -3*x + 6*y - 15, y) Out4: {y: 1, x: -3} Линейная алгебра Матрицы создаются как экземпляры специального класса: >>> from sympy import Matrix >>> Matrix([1,0, 0,1]) 0 1 В них можно использовать не только числа, но и символы: >>> x = Symbol('x') >>> y = Symbol('y') >>> A = Matrix([1,x, y,1]) >>> A x 1 >>> A**2 + x*y, 2*x [ 2*y, 1 + x*y] Над матрицами можно производить практически любые действия линейной алгебры, о чем можно ознакомиться из соответствующего раздела документации и действий над матрицами с помощью модуля mpmath. Сопоставление с образцом Сопоставление с образцом (англ. Pattern matching) — метод анализа списков или других структур данных на наличие в них заданных образцов. для проведения сопоставления используется метод .match(), с использованием объекта класса Wild. Метод возвращает словарь с соответствующим решением: >>> from sympy import * >>> x = Symbol('x') >>> p = Wild('p') >>> (5*x**2).match(p*x**2) {p_: 5} >>> q = Wild('q') >>> (x**2).match(p*x**q) {p_: 1, q_: 2} Если сопоставление не дает результатов, возвращается объект None: >>> x = Symbol('x') >>> p = Wild('p', exclude=1,x) >>> print (x+1).match(x+p) # 1 is excluded None >>> print (x+1).match(p+1) # x is excluded None >>> print (x+1).match(x+2+p) # -1 is not excluded {p_: -1} Печать Существует несколько способов того как выражения могут быть выведены на экран. Стандартный print Привычный вывод str(expression): >>> from sympy import Integral >>> from sympy.abc import x >>> print x**2 x**2 >>> print 1/x 1/x >>> print Integral(x**2, x) Integral(x**2, x) Pretty printing Формат, по которому работает isympy по умолчанию в виде псевдографического ascii-форматирования. Немного больше примеров здесь. >>> from sympy import Integral, pprint >>> from sympy.abc import x >>> pprint(x**2) 2 x >>> pprint(1/x) 1 - x >>> pprint(Integral(x**2, x)) / | | 2 | x dx | / Для того, чтобы сделать pprint по умолчанию в стандартном интерпретаторе, производим следующую процедуру: >>> from sympy import * >>> import sys >>> sys.displayhook = pprint >>> var("x") x >>> x**3/3 3 x -- 3 >>> Integral(x**2, x) / | | 2 | x dx | / Python printing >>> from sympy.printing.python import python >>> from sympy import Integral >>> from sympy.abc import x >>> print python(x**2) x = Symbol('x') e = x**2 >>> print python(1/x) x = Symbol('x') e = 1/x >>> print python(Integral(x**2, x)) x = Symbol('x') e = Integral(x**2, x) LaTeX printing LaTeX — наиболее популярный набор макрорасширений системы компьютерной вёрстки TeX, который облегчает набор сложных документов. >>> from sympy import Integral, latex >>> from sympy.abc import x >>> latex(x**2) x^{2} >>> latex(x**2, mode='inline') $x^{2}$ >>> latex(x**2, mode='equation') \begin{equation}x^{2}\end{equation} >>> latex(x**2, mode='equation*') \begin{equation*}x^{2}\end{equation*} >>> latex(1/x) \frac{1}{x} >>> latex(Integral(x**2, x)) \int x^{2}\,dx MathML MathML (от англ. Mathematical Markup Language, язык математической разметки) — это приложение XML, используемое для представления математических символов и формул в документах WWW. MathML рекомендован математической группой W3C. >>> from sympy.printing.mathml import mathml >>> from sympy import Integral, latex >>> from sympy.abc import x >>> print mathml(x**2) x2 >>> print mathml(1/x) x-1 sympy.printing Кроме того, как видно из рассмотрения python printing доступен модуль sympy.printing. Интересные функции, доступные в этом модуле: * pretty(expr), pretty_print(expr), pprint(expr): возвращают expr в pretty-формате, как показано выше; * latex(expr), print_latex(expr): возвращают expr в LaTeX представлении; * mathml(expr), print_mathml(expr): то же для MathML; * print_gtk(expr): выводит expr в Gtkmathview — GTK виджет, показывающий MathML код. Что дальше Этот урок является частичным переводом пособия для новичков. Если вам интересно, что еще может SymPy, то вы можете приступить к более полному руководству пользователя или пробежаться по модулям. Ссылки *Официальный сайт SymPy — отсюда можно скачать пакет для установки *Tutorial — основа этого урока *SymPy User’s Guide *SymPy Modules Reference *Собственное wiki Категория:Программирование и научные вычисления на языке Python